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ABSTRACT 


A  generalized  discontinuous  solution  is  found  for  the 
adiabatic  two-fluid  equations  in  the  steady  state:  it  covers 
the  case  of  strong  shocks  and  enables  one  to  give  a  complete 
account  of  the  steady  state  solutions  of  these  equations.   By 
considering  a  piston  problem  using  numerical  methods,  time 
dependent  solutions  of  the  equations  are  also  found;  these 
rapidly  steepen  and  converge  to  the  discontinuous  steady 
state  solutions  whenever  these  exist. 


2  - 


Table  of  Contents 

page 
Abstract  2 

Section 

1.  Introduction  and  Summary  4 

2.  The  Equations  of  Motion  7 

3.  Steady  State  Solutions  of  the  Equations    13 

4.  Numerical  Solution  of  the  Time  Dependent 
Equations  21 

5.  Discussion  of  the  Results  27 
Appendix  A:  Finite  Difference  Equations  28 
Table  1  :  Summary  of  the  Cases  Computed  32 
References  33 
Figures : 

1.  Domains  of  Existence  of  the  Steady  State 
Solutions  and  the  Positions  of  Cases  I 
to'T  35 

2.  Behaviour  of  the  Steady  State  Solutions 

in  the  (B,v)  plane  36 

3.  Variation  of  the  Value  of  A   just 

Behind  the  Wave  Front  in  Case  II         37 

4.  The  Wave-Form  for  X  at  t  =  60.0 
and  120.0  for  case  II;  x  is  the 
Lagrangian  Co-ordinate  38 

5.  Comparison  of  the  Time  Dependent 
Solution  and  a  Composite  Steady  State 
Solution  for  Case  III  39 

6.  Comparison  of  the  Time  Dependent 
Solution  and  a  Composite  Steady  State 
Solution  for  Case  IV  40 

7.  Comparison  of  the  Time  Dependent  and 
Steady  State  Solutions  for  CaseTT        4l 

-  3  - 


NYO-9763 


The  Development  of  Compression  Waves  In  an  Adlabatic 
Two-Fluid  Model  of  a  Collision-Free  Plasma 


1.    Introduction  and  Summary 

The  two-fluid  model  of  a  collision-free  plasma  has  received 
considerable  attention,  especially  the  case  of  zero  temperature 
■where  it  coincides  with  the  particle  picture.   Some  steady  state 
solutions  are  well-known^  ^^     '    and  Gardner  and  Morikawa^    have 
obtained  an  asymptotic  solution  for  small  amplitude  waves.   The 
simplest  extension  of  the  model  to  warm  plasmas  is  the  intro- 
duction of  an  adiabatic  pressure  and  it  has  been  suggested  that 
this  might  provide  a  first  approximation  to  a  description  of  a 
collision-free  shock  (see  Gardner,  et  al.^  '    and  Morawetz^ -^  0  • 
Banos  and  Vernon^  '   and  Vernon^  ' '   have  studied  solitary  wave 
and  periodic  steady  state  solutions  of  these  equations  and 
Gardner  and  Morikawa ' s  asymptotic  treatment  can  similarly  be 
extended  to  this  case.   Our  chief  concern  in  this  report  is  to 
study  the  development  of  a  large  amplitude  compression  wave  by 
using  numerical  methods  and  to  relate  the  situations  holding 
after  long  intervals  of  time  with  the  steady  state  solutions. 
Similar  work  using  a  different  approach  was  briefly  reported  by 
Kilb,  Auer  and  Hurwitz  in  1960^  '^^  '   and  is  now  fully  described 
in  references  ( 13 )  and  (l4). 


■X- 

This  work  was  carried  at  the  Institute  of  Mathematical  Sciences 
while  the  author  was  on  sabbatical  leave  from  the  United  King- 
dom Atomic  Energy  Authority,  Harwell. 
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The  specific  problem  that  is  studied  here  Is  the  following: 
we  have  a  half  space  of  plasma  at  rest  and  with  uniform  density, 
bounded  by  a  plane  perfectly  conducting  piston,  and  containing  a 
constant  magnetic  field  parallel  to  the  plane  of  the  piston. 
The  piston  Is  now  accelerated  Into  the  plasma  to  a  final  velocity 
at  which  It  Is  held  hereafter.   The  problem  Is  to  find  the 
ensuing  motion  of  the  plasma  and  the  changes  In  the  magnetic 
field.   Initially  we  asked  the  particular  question  -  does  the 
wave  always  break  as  In  conventional  non-disslpatlve  fluid 
dynamics?  Having  found  that  in  many  cases  this  does  occur,  we 
have  continuled  the  solution  beyond  this  point  using  the  von 
Neumann-Rlchtmyer  artificial  viscosity  to  integrate  through  the 
discontinuity  that  develops. 

The  first  result,  provoked  by  the  behavior  of  the  numerical 
solutions,  is  that  "for  sufficiently  strong  shocks"  there  exists 
a  generalized  solution  of  the  steady  state  equations  that  is 
discontinuous  in  the  fluid  variable.   This  is  described  in 
Section  3  which  contains  a  complete  account  of  the  steady  state 
solutions.   Briefly,  there  are  always  two  constant  states  and 
the  types  of  solution  that  exist  depend  on  the  compression 
ratio  r\     between  the  states  and  the  ratio   p   of  the  fluid 
to  magnetic  pressure  in  the  initial  state.   The  [t],^    )    plane  is 
divided  into  three  main  regions  --(i)  for  small  r\      the  solitary 
wave  and,  except  for  very  small  p   ,  all  the  associated  periodic 
waves  exist;  (ii)  the  solitary  wave  ceases  to  exist  when,  at 
its  peak,  the  flow  speed  exceeds   the  characteristic  speed  and 
then  only  some  of  the  periodic  waves  exist;  and  (ill)  when 
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the  characteristic  speed  exceeds  the  flow  speed  In  the  second 
constant  state,  the  discontinuous  solution  appears.   The  two 
constant  states  have  different  entropies  and  only  in  case  (ill) 
can  they  be  joined  by  a  solution. 

The  numerical  solutions  of  the  time  dependent  equations 
show  a  very  close  relationship  to  the  steady  state  solutions: 
in  case  (ill)  the  convergence  to  the  steady  state  solution  is 
very  rapid;  and  in  case  (ii),  although  convergence  is  not  so 
rapid,  a  very  close  identification  can  be  made  between  the 
actual  solution  and  one  composed  from  the  steady  state  solutions, 
This  solution  consists  of  a  partly  discontinuous  leading  edge 
followed  by  a  wave  train,  the  first  waves  of  which  are  the 
largest  periodic  solutions  that  exist  in  the  particular  case. 
In  case  (i)  the  leading  edge  appears  to  be  no  longer  discon- 
tinuous and  no  such  identification  has  been  found  possible. 

The  author  wishes  to  thank  C.  S.  Morawetz  and  J.  Berkowltz 
for  suggesting  the  problem  and  for  many  valuable  discussions. 
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2.  The  Equations  of  Motion 

We  will  take  the  piston  to  lie  In  the   (y,z)   plane  and 
let  It  "be  moved  in  the  positive  x-dlrection;  then  all 
quantities  depend  only  on  x  and  the  time   t  .   We  shall 
also  confine  ourselves  to  the  case  where  the 'magnetic  field 
remains  in  the  z-dlrectlon  .   The  only  electromagnetic 
quantities  that  will  therefore  arise  in  our  equations  are: 

B   the  z-component  of  the  magnetic  field; 

E   the  y-component  of  the  electric  field; 
and   J   the  y-component  of  the  current. 

Since  there  are  no  currents  or  electric  fields  initially, 
all  other  components  are  zero  except  for  the  x-component  of 
of  the  electric  field,  the  charge  separation  field,  which  will 
not  enter. 

Under  these  conditions  the  motion  of  the  ions  and 
electrons  Is  essentially  two-dimensional  and  the  fluid 
equations  can  easily  be  derived  from  the  collision-free 
Boltzmann  equations  for  the  two  particles  by  taking  moments 
of  their  distribution  functions  in  the  usual  way.   Further 
assumptions  that  are  made  in  this  derivation  are  that: 

(l)    there  is  charge  neutrality  --  except  in  Polss'on's 
equation  which  is  not  needed; 
and   (li)    the  pressure  tensors  of  the  ions  and  electrons 


Recent  work  by  C  S.  Gardner  has  shown  that  there  is  a  much 
wider  class  of  solutions  when  this  restriction  is  relaxed. 
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are  diagonal  and  together  yield  a  total  pressure  that  is 
isotropic  in  the   (x,y)   plane  and  satisfies  a  perfect  gas  law 
with  adiabatic  exponent  7=2. 

Thus  we  obtain  the  following  conservation  equations  where 
we  use  subscripts  to  denote  partial  differentiation: 

(2.1)     p^  +  (pu)^  =  0 


(2.2)      p(u^  +  ^  "x^  +  Px  "  "^^ 


(2.3)  e^  +  [u(e  +  p)]^  =  JE 

where 

p  =  p   +  p   ,   the  total  mass  density  of  ions  and 
electrons; 

u  =  u   =  u   ,   the  common  x-component  of  velocity; 
and   e  ,      the  total  internal  energy  density, may  be  written  in 
the  form 

(2.4)  e  =  p  +  ^  pu2  +  1   /    +  -  \  J 


2       2   [   e^e_  )  p 


m   ,  m   and   e ,  ,  e   being  the  ion  and  electron  masses  and 
charges . 

From  the  equations  for  the  conservation  of  momentum  in 
the   y-direction  we  obtain  an  equation  for  the  current: 


(2.5)     \   +  (uJ)  =  /-  !+!:i\  p(E-uB) 

y  m^m_  y 
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And  finally  we  have  two  of  Maxwell's  equations  (we  use 
rationalized,  m.k.s.  units  and  neglect  the  displacement 
current ) : 

(2.6)     B  +  E  -  0  and  B  +  (j.J  =  0 

where  [i     is  the  magnetic  permeability  of  space. 

The  initial  constant  state  is  described  by  the  three 
quantities   p  ,  p  ,  and  B   ,   the  initial  values  of   p,  p  and 
B;      and  three  functions  of  these  quantities  are  of  particular 
significance : 

B 

a  =  T-77:  ,  the  initial  Alfve'n  speed: 

o    /    \l/2  '  I-    J 

P   =  — ;r ,   the  ratio  of  the  initial  fluid  and 

°'  (B^/2^L) 

magnetic  pressures; 

/  e^e_    \l/2         -1/2 

and  X   =  (  -  M-P  1   =  ( (i-k  co  )      ,      which  can  be 

o    I    ni,^_    oj  p'      ' 

-I/2 
regarded  as  either  the  ratio  of  the  speed  of  light   (m-  k) 

to  the  plasma  frequency  cu  ,   or  as  the  geometric  mean  of  the 

gyromagnetic  radii  of  ions  and  electrons  moving  across  the 

initial  field  with  speed  a 

^      o 

■X- 

Introducing  first  the  difference  velocity 


In  fact,  the  difference  between  the  ion  and  electron 
velocities  is^       \       / 

r.  /   e  ,e   \  1/  2  /  m^   m 
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(2.7) 


mm   \  2 


^+^. 


J 

P 


In  place  of  the  current,  we  normalize  all  our  variables  In 
terms  of  these  fundamental  quantities  in  the  following  way; 


X 

u 

p 

p 


o 


BiM 


t 

X   /a 

0^      0 

V    = 

=   a 
o 

B 

.    B 

0 

e 

B^/. 

E  :  a  B 
o  o 


There  will  be  no  confusion  if  we  now  denote  the  normalized 
variables  by  the  same  symbols  as  before.   Equation  (2.5) 
becomes 


(2.8) 


V,  +  u  V   =  E-uB 
t       X 


which  on  combination  with  (2.6)  yields  an  equation  for  B  +  v 
identical  with  the  continuity  equation  for  p  :   since  these 
quantities  are  equal  initially,  they  remain  so.   In  addition 
we  may  substitute  for  the  internal  energy 


(2.9)     e   =  ^   p(u^  +  v^)  +  p 


in  the  energy  conservation  equation,  obtaining  the  usual 
relation  that  the  entropy  is  constant  along  the  particle  paths. 
Thus  our  equations  become: 
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(2.10)     p^  +  (up)^  -  0 


(2.11)     u^  +  uu^  +  (l/p)p^  -  vB 


(2.12)     A,  +  uA   =  0 


(2.13)     v^  =  p-B 


(2.14)     pv  +  B^  -  0  ; 


where 


(2.15)         P  =  Ap 


2 


The  initial  conditions  are  now 


p=B-l     ,      A-p  =  |p^ 
(2.16) 

u  =  V  =  0 


The  piston,  whose  velocity  will  be  prescribed,  we  shall  assume 
is  a  perfect  conductor.   Thus   E  will  equal   uB   just 
inside  the  plasma  and,  from  (2.8),  the  derivative  of   v   along 
the  piston  path  will  be  zeroj  our  boundary  conditions  on  the 
piston  become 


(2.17)    u  =  Up(t),  v  =  0 
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It  will  be  noticed  that  when  p   =0  these  boundary  conditions 
together  with  equation  (2.II)  will  not  allow  the  piston  to  be 
accelerated.   For  p  >  0  a  boundary  layer  will  develop  at  the 
piston  which  will  become  Increasingly  severe  as   p  — ^  0  and 
prevent  us  from  computing  the  limiting  case. 

Equations  (2.IO)  to  (2.I5)  a^^e  mixed  hyperbolic-parabolic, 
their  characteristics  being  the  usual  characteristics  forn  non- 
dlsslpatlve  fluid  flow  plus  two  parallel  to  the  x-axls  which 
are  the  limits  of  the  characteristics  of  the  electromagnetic 
wave  equation  as  the  speed  of  light  tends  to  Infinity.   When 
the  plasma  Is  assumed  perfectly  conducting,  by  (2.8)  the 
derivatives  of  v   can  be  neglected,  so  that   p  =  B   and  the 

equations  reduce  to  those  of  ordinary  fluid  dynamics  with  the 

1   2 
"total"  pressure   p  +  -p  B   . 

Because  of  the  relevance  of  this  "one  fluid"  theory  In  at 

least  parts  of  the  flow  one  might  expect  compressive  waves  to 

steepen  and  possibly  form  discontinuities.   The  form  of  these 

Is  apparent  from  the  equations:   If   p   and  B  remain  finite, 

then  by  (2.I7)  so  will   v  ;  thus  there  can  be  no  jump  In  B 

or  In  V   at  the  discontinuity,  which  will  Involve  only  the 

fluid  quantities   p,  u,  and  p  .   The  situation  Is  then  different 

from  that  In  the  one-fluid  theory  so  perhaps  one  should  not 

always  expect  continued  steepening  to  the  point  of  breaking. 

The  Jump  conditions  at  any  discontinuity  that  does  arise, 

however,  will  be  just  the  Ranklne-Hugonlot  relations  In  the 

usual  fluid  variables.   This  does  not  mean  of  course  that  the 
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shock  that  this  might  correspond  to  In  a  more  accurate  model 
would  necessarily  be  an  ordinary  viscous  shock. 


3.    Steady  State  Solutions  of  the  Equations 

In  this  section  we  shall  consider  solutions  of  the 
equations  (2. 10)  to  (2.I5)  which  are  steady  when  viewed  from 
a  frame  of  reference  moving  with  a  constant  velocity   U  ,  the 
"shock"  velocity.   The  existence  of  solitary  wave  and  periodic 
solutions  has  been  demonstrated  previously  for  both  cold  and 
warm  plasmas^  '  '  '       .      However,  we  cannot  expect  these  to 
completely  describe  the  situation  in  our  problem  long  after 
the  piston  has  attained  its  final  velocity,  for  they  do  not 
satisfy  all  the  boundary  conditions;  Indeed,  there  cannot  be 
any  continuous  steady  state  solution  of  these  non-dlsslpatlve 
equations  joining  initial  and  final  states  with  differing 
entropies.   Thus  we  shall  need  to  look  for  both  continuous 
and  discontinuous  solutions;  we  shall  see  later  that  both  play 
roles  in  the  behavior  of  our  system  so  we  shall  try  to  give  a 
complete  description  of  all  solutions  in  this  section. 

Attaching  primes  to  quantities  referred  to  the  moving 
frame  of  reference,  the  conservation  equations  become  on 
Integration 

(3.1)     pu'  =  -U 

(3-2)      -Uu'  +  p  +  1/2  B^  =  U^  +  1/2  p   +  1/2 
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and 


(3.3)      -u'(e'  +  p)  +  UB  =  Up^  +  I  U^  +  U 


where  the  constants  of  integration  have  "been  determined  from 
the  Initial  state;  in  the  last  equation  we  have  used  the 
deduction  from  (2.6)  that   E  =  E  -  UB  =  constant  =  -  U.  .  Now  we 
replace  the  density  p   by  the  specific  volume   V  =  —  and 
substitute  into  (3o)  the  expressions  for  u   obtained  from 
(3-1);.  for   p   from  (3.2),  and  for   e   from  (2.9);  the  result 
is  a  quadratic  for  V  in  terms  of  B  and  v  : 


;3.^)    3U^V^  -  2V[2U^+  1  +  P  -  B^]  +  [U^-  v^-  2B  +  2  +  2p  ]  =  0 


The  remaining  equations  (2.1;?)  and  (2.l4)  become 


(3.5)  V  g  +  BV  -  1  =  0 


and 


dB 


(3.6)  V  ^  +  V  =  0 

dx 


Thus  solutions  are  described  by  this  pair  of  equations 
for  B  and   v  with  V  defined  by  one  of  the  roots  of  (3.4). 
Constant  states  correspond  to  singularities  of  the  differential 
equations,  given  by  v.  =  0  ,   B.V.  =  1  .   If  we  substitute 
1/B  for  V  in  (3.4)  with  v  =  0  ,   we  find  one  root 
corresponding  to  the  initial  state   B  =  V  =  1  ,   and  the 
other,   B-.  =  ri  =  l/"^-!  j   which  we  would  like  to  identify  with 
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a  final  state  and  for  which 
(3.7)     u2  = 


3-11 


Clearly,  when  the  compression  ratio  ri   is  greater  than 

unity,  the  shock  velocity  U   is  greater  than  the  magneto-sonic 

speed   (1  +  p^)^  ^  . 

The  solutions  are  most  conveniently  described  in  the 
(B,v)  plane;  with  these  co-ordinates,  equations  (3-5)  and  (3-6) 
define  a  direction  field  on  a  two-sheeted  surface,  each  sheet 
corresponding  to  a  root  of  the  quadratic  {^.h)    for  V  .   The  two 
sheets  are  joined  on  a  closed  convex  curve  on  which  these  two 
roots  are  equal.   Inside  the  curve  the  roots  become  complex  so 
that  no  physically  admissible  solution  may  cross  into  this 
region.   If  we  write  the  equation  for  the  curve  in  the  form 

\  „   3(1+P^)(n-1)(  2    12(l+p^)B(B-r|) 

(3.8)    |b2-_-3-^ -J   -  3-^ -^3UV^0 

we  see  that  it  always  lies  to  the  right  of  the  second 
singularity   (B  ^   r\,    v     =   O);      it  passes  through  this  point 
for  the  particular  value  of   -q   given  by 


:3.9)  -|ii.i^  =  -i 


X\      (3-Tl)       1+pQ 

Two  other  physical  conditions  may  further  confine  the  position 
of  solutions  on  the  two  sheets:   firstly,  the  entropy  cannot 
fall  below  its  initial  value,  i.e.,   A  >  p^/2  ;      and  secondly 
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the  density  must  remain  positive,  i.e.,  V>  0.  However, 
when  p  is  replaced  by  a/v  the  momentum  equation  (3. 2) 
can  te  written  as  a  cubic  for  V  , 

(3.10)  ■    V^  +  V^  ).  ^  "  ^  "  ^o  -  1  V  +^=   0 

the  two  positive  (or  complex)  roots  of  this  being  those  of 
{^.^).  Thus   only  when  A  =  p   =0   does   V  become  zero  and 

then  the  two  roots  become  zero  together  so  that  the  condition 
that  we  remain  outside  (3-8)  already  contains  this  restriction. 

Now  let  us  turn  to  the  nature  and  position  of  the 
singularities.   The  discriminant  for  the  characteristic  equation 
takes  on  the  simple  form  V~  5(VB)/^B  at  the  singularities, 
which  on  computing  the  derivative  becomes 

M?-l-p. 

(3.11)  A   =   \  g  \    ,  (1  =  0,  1)  . 

^    V  (M^-p  ) 


Here 


M?  =  u.'^p.  B?  =  V-?U^,   the  square  of  the  Alfven  Mach 


number,  and 

p 
6.  =  2p.  Bt  ,   the  ratio  of  the  fluid  to  the  magnetic 

p    '^ 
pressure.   At  the  initial  state  M  =  U"  >  1  +  P^   for  n  >  1 

so  that  this  singularity  is  always  a  saddle  point  in  the  cases 

we  are  Interested  in.   It  is  also  easily  ascertained  that  it 

always  lies  on  the  sheet  corresponding  to  the  larger  of  the 

2 
roots  for  V  .   At  the  second  singularity  M   is  always  less 

2 
than   1  +  P-,   so  that  the  singularity  is  a  center  if  M-j^  >  (3^ 
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2 
and  a  saddle  point  If  M  <  p   .   This  condition  can  be  written  in 

the  form 


(5.12)  ^^     >  ^ 


2/^   ^   1+P 

If   r|   is  so  small  that  the  condition  is  not  satisfied  then  the 
singularity  is  a  center  and  also  it  lies  on  the  sheet  of  the  larger 
root;  if  (3,12)  is  satisfied  then  it  is  a  saddle  point  and  it  lies 
on  the  sheet  of  the  smaller  root.   The  changeover  occurs  when  r| 
satisfies  (3-9)  when,  as  we  have  already  seen,  both  sheets  contain 
the  singularity.   Thus  a  solution  can  always  leave  the  initial 
state  on  the  sheet  of  the  larger  root  but  only  when  (3. 12)  is 
satisfied  can  it  enter  a  final  state  and  this  is  on  the  other 
sheet.   To  make  the  change  from  one  sheet  to  the  other  there  must 
be  a  discontinuity  in  V  ,   which  will  occur  where  the  two 
solutions  issuing  from  the  end  states  meet  in  the   (B,v)  plane. 
This  is  what  was  expected  from  the  time  dependent  equations.   It 
is  worth  noting  that  to  obtain  a  change  from  one  sheet  to  the 
other  which  is  continuous  in  V  would  involve  both  solutions 
being  tangential  to  the  curve  (3-8);  this  does  not  occur  in  general, 
but  one  particular  case  is  of  interest  as  we  shall  see  in  the  next 
section . 

We  are  now  in  a  position  to  tie  all  these  observations 
together  in  a  complete  description  of  the  three  main  types  of 
steady  state  solution.   For  any  given  value  of  p  (O  <  p   <  0°) 
the  situation  depends  on  one  other  parameter,  either   U  or  t] 
which  are  connected  by  equation  (3-7)j  we  shall  use  the  latter 

~  17  - 


which  lies  in  the  range   ( 1  <  t]  <  3 )  j  and  the  domains  in  which 
each  situation  holds  are  shown  in  Figure  1.   For  sufficiently 
weak  disturbances,  i.e.,  x]      near  unity,  a  solitary  wave 
solution  exists  and  appears  as  a  lobe  in  the   (B,v)   plane 
about  the  center  at   B  =  t]  ,   v   =  0  j   inside  this  is  a 
complete  set  of  periodic  solutions  of  increasing  entropy  and 
decreasing  amplitude  as  they  circle  nearer  and  nearer  to  the 
center.   This  situatior)  (which  is  illustrated  by  two  of  the  cases 
shown  in  Figure  2)   holds  because  either  the  domain  D,.  ,  in 
which  the  roots  of  {^.^)   are  complex,  is  empty  or  it  lies  to 
the  right  of  the  solitary  wave  lobe.   D-^  is  empty  in  the  region 
above  and  to  the  left  of  curve   C-,   in  Figure  1;  on  this  curve 
it  appears  as  a  single  point  which  for  p  >  0.03704  (or 
T)  <  1.6)   is  outside  the  solitary  wave,  while  otherwise  it  is 
inside.   In  this  latter  case,  as   t]   is  increased   D   expands 
excluding  more  and  more  of  the  periodic  solutions  until  it 
intersects  the  solitary  wave  lobe;  this  region  is  shown  as  a 
small  triangle  bounded  by  curves   C  ,  C   and  the  line  p   =0 
in  Figure  1  .   For  the  larges  values  of  p    the  domain  D„ 
approaches  the  solitary  wave  lobe  from  the  outside;   in  both 
cases  Intersection  occurs  on  the  curve   C   given  by 


(3.13)  ^Li^-.l 


o 


and  beyond  this  the  solitary  wave  no  longer  exists.   Both  the 
solitary  wave  and  the  periodic  solutions  lie  entirely  on  the 
sheet  corresponding  to  the  larger  root  for  V.   The  other  sheet 
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contains  no  singularities  and  Its  solutions  consist  of  a  series 
of  curves,  each  crossing  the   B-axls  at  just  one  point  and 
curving  to  the  left. 

As  T\      Increases  beyond  the  value  given  "by  (3-13),  the 
domain  D„   approaches  the  second  singularity  (B  ,,  v.,), 
excluding  more  and  more  periodic  solutions,  until  on  curve   C^ 
Its  boundary  passes  through  It.   Beyond  this  point   ( B-,  ,  v,  )   is 
a  saddle  point  on  the  second  sheet.   In  this  case.  Illustrated  by 
the  last  diagram  of  Figure  2,  there  are  two  particularly 
interesting  solution  curves,  one  issuing  from  the  first  singu- 
larity on  one  sheet  with  positive   dv/dB,   and  one  from  the 
second  on  the  other  sheet  with  a  negative   dv  dB:   they  meet  at 
a  point  on  the   (B,v)   plane  where  a  jump  discontinuity  in   V 
joins  them.   Plotted  against  the  Lagranglan  space  co-ordinate 
they  will  appear  as  in  Figure  J.      This  is  the  only  steady  state 
solution  that  joins  two  distinct  constant  states. 

The  above  description  has  been  almost  entirely  mathematical 
but  the  physical  meaning  of  the  restrictions  placed  on  each 
type  of  solution  can  now  be  made  clear.   The  gas  sound  speed, 
which  is  also  the  characteristic  speed  is   (2pV)  ^         and  in  the 
constant  states  equals   (P.B.)  /'  ;  also  the  flow  speed. is 

u'.    =   M.Bt/^  .   Thus  the  condition  M-,  <  B,   for  the  discontinuous 
111  1  ^  ^1 

solution  to  exist  is  merely  the  usual  shock  condition,  that  the 
flow  behind  the  shock  be  subsonic,  applied  at  the  constant 
state.   That  there  is  such  a  gap  before  this  solution  "exists  is 
partly  due  to  the  role  played  by  the  one-fluid  model,  in  which 
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l/? 

the  characteristic  speed  is  [(l+3^)B.]  '  .  In  the  two  constant 
states,  BV  =  1  and  v  =  0  and  the  conservation  equations  (3.1) 
to  (3-3)  applied  to  these  states  reduce  to  the  Rankine-Hugoniot 

conditions  for  the  one  fluid  model  with  a  total  pressure 

12  -    1   2 

p  +  :5-  B   and  total  energy  e  +  —  B  V.   In  particular,  equati 


p  ^  l^llVt  WWWl^_L  ^J.1V_J.     gjjr  ^  .  p 


ion 


(3-7)  for  the  shock  speed  corresponds  to  Prandt ' s  shock  relation 


(3.1^)    H^u^^  +  (l-H^)c^  =  u^u^ 


2    1        2 
where   u  =  ^  and   c   =  l+p   .   Thus  one  arrives  at  the  twin 
3        o      o 

conditions  for  the  existence  of  the  discontinuous  solution, 

2  2 

M  >  1+P    and  M  <  p   ,   which  clearly  leaves  a  gap  in  the 

range  of  values  of  t\      for  which  they  can  be  satisfied. 

In  addition,  Vernon^  '  ^  expresses  the  condition  (3,13)  for 

the  limit  of  the  domain  in  which  the  solitary  wave  exists,  as 


(3.15)         V  =  V     where   v-^  =  p  /U^ 


and  V   is  the  value  of  V  at  the  peak  of  the  solitary  wave. 

This  condition  can  be  written  in  the  form  u  =  VU  =  (f3  /V) 

1/2 
=  (2pV)     ,  i.e.,  that  for  an  infinitely  weak  shock.   We  shall 

see  in  the  next  section  that  in  the  domain  between  curves   C 

and   C    of  Figure  1   discontinuities  do  occur  even  though 

they  cannot  directly  link  constant  states. 
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^.    Numerical  Solution  of  the  Time  Dependent  Equations 

For  the  numerical  calculations  It  Is  most  convenient  to 
work  with  the  equations  In  their  Lagranglan  form.   In  addition, 
to  enable  the  computations  to  be  simply  continued  through  the 
shock  fronts  we  Introduce  a  von  Neumann-Rlchtmyei-  ^ -^"^^  artificial 
viscosity.   The  equations  are  then 

(^.1)  u,  +  [^  +  q  +(^)b2]  =  0 


(4.2)  V,  =  u 

t     X 


(4.3)  A^  +  (qV)V^  =  0   where 


I  qV  =  a^h^(u  )^   If   u  <  0 


I qV  =  0  If   u   >  0 

'v. 


(4.4)  B    =  BV  -  1 


The  parameter  a,  which  Is  normally  unity,  enables  us  to  make  a 

simple  transformation  of  the  variables  to  deal  with  cases  of 

2 
very  high  (3  ;   and   a   Is  a  parameter  used  to  control  the 

amount  of  artificial  viscosity  Introduced.   The  Initial  and 

boundary  conditions  are, 

(4.5)      t  =  0,   X  >  0:   B  =  V  =  1,   u  =  q  =  0,   A  =  p  /2 


O' 


(4.6)      t  >  0,   X  -  0:   B   =0,   u  =  Up(t 
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The  difference  scheme  used  to  approximate  these  equations 
is  fully  described  in  Appendix  A;  it  allows  explicit  integration 
of  the  first  three  equations  while  {h.k)    is  solved  implicitly; 
all  differences  are  centered  on  a  staggered  mesh  except  for  the 
expression  for  q   in  the  first  equation.   Criteria  for  the 
stability  of  the  difference  scheme  are  effectively  the  same  as 
those  for  pure  fluid  dynamics, there  being  one  for  the  shocked 
region  and  one  for  the  smooth  flow  outside  this  as  described, 
for  example,  by  Richtmyer    .  A  running  check  on  the  computation 
is  provided  by  the  energy  conservation  integral, 

CO 

(4.7)   Up[p  +  q  +(|)b2  ]p=  ft  l"  fi  ^^+  P^)  +  "^i  ^x+  2  ^^V)]d^ 

where  the  suffix   P   is  used  to  refer  to  qualities  evaluated 
at  the  piston. 

For  Illustrative  purposes  attention  has  been  concentrated 
on  five  cases  for  each  of  which  p   =  0.1;  in  each  case  the 
piston  is  accelerated  uniformly  from  rest  to  a  final  velocity 
in  either  one  or  two  units  of  time;  it  is  then  held  at  this 
velocity  thereafter,  the  five  cases  corresponding  to  final 
velocities  of  0.1,  0.4,  0.5,  1.0  and  1.5.   After  an  initial 
period    in  which  the  effect  of  the  acceleration  is  dominant, 
the  solutions  in  general  settle  down  to  a 'fairly  common 'form 
which  consists  of  four  portions.   The  disturbance  propagates 
into  the  plasma  at  a  more  or  less  definite  speed  and  with  a 
leading  edge  that  may  steepen  and  break  into  a  partly 
discontinuous  form.   This  is  followed  by  a  wave  train  of 

-  22  - 


increasing  length,  the  amplitude  of  the  waves  decreasing  away 
from  the  disturbance  front  until  a  constant  state  is  reached; 
the  constant  state  makes  up  the  third  portion  of  the  solution 
and  extends  almost  to  the  piston  to  which  it  is  joined  by  a 
narrow  boundary  layer.   Figure  4  shows  a  typical  example. 

For  at  least  three  of  the  cases  a  detailed  description 
of  the  results  can  be  obtained  by  reference  to  the  steady  state 
solutions  described  in  the  previous  section.   The  only  difficulty 
arises  in  the  identification  of  each  case  with  a  point  in  Figure 
1,  which  determines  what  steady  state  solutions  are  available. 
The  matching  of  p    is  immediate  and  for  the  other  parameter 
we  have  identified  the  final  piston  velocity  with  the  fluid 
velocity  in  the  second  constant  state  and  used  the  steady  state 
equations  to  obtain  r|,U  etc.   Another  possibility  would  be  to 
measure  the  shock  velocity  U  from  the  numerical  results  and 
use  this  as  the  starting  point  of  the  natching;  but  we  prefer 
to  attempt  a  complete  prediction  from  given  data  even  if  this 
does  have  a  greater  error.   Using  the  first  method,  then,  the 
five  cases  correspond  to  the  points  m.arked  with  a  cross  in 
Figure  1  and  the  situation  on  the  (B,v)  plane  for  each  case  is 
shown  in  Figure  2. 

As  would  be  expected,  case  "Y^    is  the  most  straightforward; 
a  discontinuous  steady  state  solution  joining  the  initial  state 
to  the  state  corresponding  to  the  piston  velocity  exists  and 
our  matching  is  therefore  correct.   The  piston  reached  its 
final  velocity  at   t  =  2,0   and  by  t  =  12-5  the  front  of  the 
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wave  looks  as  In  Figure  J;    the  corresponding  steady  state 
solution  Is  shown  on  the  same  figure.   There  Is  very  close 
agreement  except  In  a  region  one  unit  wide  In  which  the 
artificial  viscosity  resolves  the  discontinuity  In  V.   The 
shock  speed  and  the  state  behind  the  front  agree  very  well  too, 
as  can  be  seen  from  Table  1. 

In  case  IV  the  wave  train  behind  the  discontinuity  appears 
and  In  case  III  these  waves  are  of  Increased  amplitude  and  wave- 
length ;  parts  of  the  solutions  (at  t  =  20  and   t  =  40, 
respectively)  are  shown  In  Figures  6  and  5-   By  this  time  the 
leading  part  of  the  solution  has  reached  a  steady  form  and  the 
first  wave  of  the  following  train  Is  beginning  to  have  a  readily 
Identifiable  wavelength.   Turning  to  the  (B,v)  diagrams  In 
Figure  2,  we  can  see  how  this  form  Is  made  up  from  the  steady 
state  solutions.   The  second  singularity  Is  still  on  the  sheet 
corresponding  to  the  larger  root  for  V  and  around  It  there  Is 
a  set  of  periodic  solutions,  that  with  largest  amplitude  touching 
the  boundary  of  the  domain  D-y  at  a  point  on  the  B  axis.   Also 
passing  through  this  point  and  touching  D   Is  a  solution  on 
the  other  sheet,  using  the  smaller  root  for  V,  which  eventually 
Intersects  the  curve  corresponding  to  what  Is  left  of  the  solitary 
wave.   Where  these  two  solutions  and  the  boundary  of   D„  touch 
Is  the  only  place  where  a  solution  may  change  from  one  sheet  to 
the  other  without  causing  a  discontinuity  In  V  and   A.   Thus 
a  solution  can  be  made  up  of ( 1 )  an  Initial  part  along  the 
solitary  wave  curve,  (ll)  a  jump  to  the  second  sheet,  (ill) 
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continuation  on  this  sheet  to  the  point  of  contact  with  the 
boundary  of  D^  on  the   B-axls,  ending  with  (iv)   the  continuous 
traversal  of  the  periodic  solution  loop.   Solutions  constructed 
In  this  way  are  shown  for  cases  III  and  IV  on  Figures  5  and  6 
below  those  obtained  for  the  time  dependent  problem.   There  Is 
a  very  considerable  agreement  even  at  this  early  stage  In  the 
development  of  the  wave  train  and  with  the  matching  of  parameters 
described  above. 

In  Table  1  the  shock  velocities  and  the  values  of  the 
variables  in  the  constant  state  that  develops  near  the  piston 
are  compared  with  those  obtained  from  the  steady  state.   As  we 
should  expect,  the  shock  velocity  and  the  characteristics  of 
the  wave  train  are  consistent,  their  difference  from  the  steady 
state  values  arising  mainly  from  incorrect  choice  of  the 
parameter  x]      ;    but  the  constant  state  cannot  be  properly 
compared  with  anything  from  the  steady  state  solutions  for  its 
entropy  is  that  of  the  largest  periodic  solution  forming  the 
head  of  the  wave  train,  and  not  that  of  the  constant  state  at 
the  second  singularity.   Further,  it  is  joined  to  the  front  of 
the  disturbance  by  the  time  dependent  part  of  the  wave  train. 

The  values  of  the  parameters  for  case  II  suggest  that  a 
different  type  of  solution  may  be  obtained  in  which  there  are 
no  discontinuities  and  the  solitary  wave  plays  a  larger  part. 
In  fact,  we  have  been  unable  to  describe  the  solution  in  terms 
of  the  steady  state  solutions  in  either  this  case  or  in  case  I. 


_  p 
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The  development  of  the  solution  in  case  II  is,  however,  rather 
interesting.   This  is  illustrated  by  Figure  3   in  which  we  have 

plotted  the  value  of  A   at  the  first  minumum  in  V   (i.e., 
just  behind  what  would  be  the  shock)  as  a  function  of  time. 
Apart  from  an  early  excursion  near  the  time  that  the  piston  is 
being  accelerated,   A   increases  gradually,  due  to  steepening  of 
the  front  coupled  with  the  artificial  viscosity,  until   t''^-45; 
it  then  rises  sharply  as  the  wave  appears  to  break,  reaching  a 
maximum  at   t  =  60;  from  here  it  decreases  at  a  steady  rate  as  the 
break  weakens  until  from  t  =  105   to   t  =  I50  it  remains  at  a 
uniform  level  corresponding  to  a  very  steady  leading  edge  of  the 
wave.  The  whole  of  the  wave  train  at   t  =  60  and  the  leading 
part  at   t  =  120  are  shown  in  Figure  4;  the  most  significant 
measurable  parameter  seems  to  be  the  wavelength  of  the  leading 
waves  which  remains  at   10. 7  from   t  =  100  .   No  quantitative 
explanation  of  this  has  been  found  but  the  similarity  of  the 

problem  with  that  of  bores  suggests  that  a  treatment  like  that 

(12 
given  by  Benjamin  and  Lighthill   ^ might  be  possible. 

The  low  piston  speed  for  case  I  was  chosen  with  the  hope 

that  the  behavior  of  the  solution  might  approximate  the 

asymptotic  form  given  by  Gardner  and  Morikawa   for  a  cold  plasma. 

The  equations  that  they  give  can  also  be  derived  for  the  present 

-1  /2 
case  and  implies  that   1-B,   1-p,   and  u(l-i-p  )~     are  all 

equal;  moreover,  in  the  particular  solution  they  give  the  waves 

have  a  definite  amplitude  related  to  the  piston  velocity  and 

1/3 
widen  with  time  like   t     .   At   t  =  100  when  the  integration 

of  case  I  was  stopped  the  amplitude  of  the  waves  was  still 
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growing  although  It  was  already  much  greater  than  that  predicted, 
None  of  the  other  features  were  reproduced  to  a  useful  degree 
of  approximation. 


5.    Discussion  of  the  Results. 

We  have  shown  that.  In  the  present  model  when  the  shock 
speed  or  compression  Is  too  great  for  the  solitary  wave  to 
exist,  there  are  still  steady  state  solutions.   These  involve 
discontinuities  at  which  entropy  Is  lost  and  lead  either  to 
periodic  waves  with  a  distinctive  form  or.  In  the  case  of  the 
greatest  shock  speeds,  to  a  constant  state.   Moreover,  in 
these  cases  the  solution  of  the  time  dependent  piston  problem 
converges  to  these  solutions. 

The  mechanism  by  which  entropy  is  lost  at  the  discontinuities 
is  completely  unspecified  so  that  we  have  so  far  refrained  from 
calling  these  phenomena  colllslonless  shocks.   The  distance 
scale  on  which  the  discontinuity  appears  is  the  familiar        -.  , 
geometric  mean  of  ion  and  electron  Larmor  radii  so  its  resolution 
demands  a  mechanism  operating  on  a  scale  length  related  to 
either  the  electron  Larmor  radius  or  perhaps  the  Debye  radius. 
It  may  be  significant  that  on  calculating  the  charge  separation 
field  from  the  solitary  wave  solution,  Vernon  found  that  it 
became  Infinite  where  the  solitary  wave  ceased  to  exist.   To 
neglect  this,  as  we  have  done,  is  clearly  unjustified  and 
should  be  the  first  thing  Included  in  a  more  detailed  model  of  . 
the  shock  structure.  _„ 


Appendix  A:   Finite  Difference  Equations 

Using  a  constant  mesh  size  h  and  a  time  step  k  ,  let 
u  be  correctly  centered  at  points  x  =  jh  ,  t  =  (  n  -  -Djm 
(j,h  =  0,1,...)   and  denote  its  value  at  these  points  by  u.  j 

J 


simi 


larly  let  V,B,   and  A   be  centered  at  x  =  (j  +  ^)h  , 


n 
t  =  nk,   their  values  being  denoted  by  V.   etc.,  finally  qV 

(J 

and  q   will  be  centered  at  x  =  (j  +  ^)h,   t  =  n  -  ^)k  .   We 
then  have  the  difference  equations 


/  .    -,  \      n+1    n     k  )   J+l        i     ,   n      n 
J+l     J+l    h    /  n   \2     /vn>?    ^J+1    ^J 


-I '(^5.1) 


:a..)      v-i.v^.^    u-J-u-1 


(A. 3)  Af  1   =   A"   -    (qV)f  l(vfl    -   v"  ) 

J  J  J  J  J 


2        („.,  „,,) 


/    TrNn+l  2    \    n+1  n+ll     .„]   n+1        n+1  i     .    ^ 

where      (qV)  .        =   a      <)u.^^    -   u.      ^    if  <ju.^^-   u.      |  <   0 

„  . „  I    n+1        n+1  ,   ^    ^ 


(A.  4)  E^+j    -    (2    +  h2v^+l)B^+l    +   B^.l]    =    -Y? 
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and 

2 


(A.  5)      qr"  =    (qV)\! 


j        '^  ^j    ^  v^+^  +  v^ 


These  equations  are  solved  In  the  order  given.   When 
solving  (A.2),   u     Is  obtained  from  the  expression  for  the 
piston  velocity  which  Is  given;  and  (A. 4)  which  Is  Implicit 
needs  boundary  conditions  at  each  end.   At  the  piston  we  have 
B^-,   =  b'^    while  the  other  boundary  Is  allowed  to  move 
outwards  by  adding  points  until 


^jS  -  i-o<  ^ 


where   £   Is  some  prescribed  error.   Since  V  decays  to  Its 
value  at  x  =  +  oo  more  rapidly  that   B,   B-1.0  has  a  tall  of 
the  form  e""^  and  we  use  this  as  the  outer  boundary  condition. 
The  trl-dlagonal  matrix  equations  for  B  are  solved  In  the 
usual  recurrence  relation  manner  (see,  for  example,  Rlchtmyer  ^) 
For  the  more  accurate  equations  which  we  usually  used 

=   -  h2 


where  we  have  dropped  the  superflxes,  the  procedure  becomes: 


E^  =  1   ,        P^  =  0 


29  - 


(A. 7) 


h 


E. 


J+1 


[(2  +  5h^/6   V  )-(l  -  h2/l2   V    )E  ] 


P 


h^  +  (1 


12  M 


,]  - 1   , 


J+1 


[(2  +  5hV6  V  )-(l  -  hVl2  V    )E  ] 


j  =  1,  2,  ...,  J 


followed  by 


(A. 8) 


where 


(A. 9) 


B 


^j+i  +  n  - 1 


J+1 


n  = 


r-  E 


J+1 


1/2 


1  + 


12 


h' 


+  h(l  +^) 


B.  =   F.^^   +E.^^B.^^  ;    J.J,  J-1,  ...,  0. 


Stability.      In  regions  of  smooth  flow  away  from  shocks  we 
may  take  q  =  0  .    Then  the  amplification  matrix  for  equations 
(A.1)-(A.3)^  (a. 6)  in  terms  of  the  variables   u,  V,  and  A   is 
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2 
21d 


(A. 10)    I     21g       l-4g^d^ 


0 


where  g  =  (k/h)  sin  -^  mh  =  (k/h)s 

P     IP 
p    PA      p     h(l-^sj 

and  d^  =  ^  +  aB^  — 5 — 2 — 

V^        43^  +  h^(l  -  I  s2)V 

=  2Ayv^  +  O(h^),    normally. 

The  corresponding  characteristic  equation  Is 


(A. 11)  (A-l)[A^  -  (2-4g^d^)A  +  1]  -  0 


so  that  for  the  most  practical  purposes  the  stability  condition 
is  the  usual  one  that   (h/k)  >  {2k N    )   '     ,    the  sound  speed  (in 
Lagranglan  coordinates  ) . 

In  the  shocked  region  the  magnetic  field  and  Its  first 
derivative  are  continuous  so  that  In  this  narrow  region  the 
equations  are  Identical  to  those  of  fluid  dynamics:   for  these, 
Rlchtmyer   gives  the  condition: 


k  /2A Y  ^  [n(n  :  1/3)] 


1/2 


h  (  3J  ^  ^"^2a(Ti-l'3 


where   t]   Is  the  shock  stength,  or  compression  ratio  of  the 
shock.  _  31  _ 
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FIG.  5.    COMPARISON     OF    THE    TIME     DEPENDENT     SOLUTION 
AND   A    COMPOSITE      STEADY     STATE     SOLUTION     FOR    CASE    JR. 
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